knitr::opts_chunk$set(
message = FALSE,
warning = FALSE,
results=FALSE,
fig.align="center"
)
library(tidyverse)
library(tidycensus)
library(sf)
library(kableExtra)
library(patchwork)
library(ggplot2)
library(RSocrata)
library(viridis)
library(spatstat)
library(raster)
library(spdep)
library(FNN)
library(grid)
library(gridExtra)
library(knitr)
options(scipen=999)
options(tigris_class = "sf")
options(tigris_use_cache = TRUE)
source("functions.r")
p5 = plasma(5)
p2 = plasma(2)
crs=7131
According to FBI UCR data, San Francisco is one of the most crime-ridden city in 2019, with a extremely high property crime rate of 5506.0 crimes out of 100,000 people, which is 2.61 times the U.S. average.
This high property crime rate was making residents in San Francisco feel unsafe1 , and may prevent migrant population growth and potential population loss from dangerous areas.
read.csv("./data/crime_by_city_2019.csv")%>%
mutate(City=gsub("(Metropolitan|\\d|Police Department)","",City),
State=na_if(State,""),
Population=as.numeric(gsub(",","",Population)),
PropertyCrimeRate =as.numeric(gsub(",","",Property.crime))/Population*100000)%>%
fill(State)%>%
filter(Population>500000) %>% na.omit()%>%
mutate(City=reorder(City,PropertyCrimeRate,mean),
toHighlight=ifelse(City=="San Francisco","",0),
toHighlight=na_if(toHighlight,0))%>%
ggplot(aes(PropertyCrimeRate,City))+
geom_col(aes(fill=toHighlight),show.legend=F)+
geom_vline(aes(xintercept=2109.9),linetype=2)+
annotate("text", x=2200, y=1,label="US AVG",size=3,angle=90,vjust=1,hjust=0)+
scale_fill_manual(values=p5[3], guide="none",na.value="grey75")+
labs(x="Property crimes per 100,000 population",
y="Cities(population>500,000)",
title="Property Crime Rate of US Major Cities",
subtitle="FBI UCR 2019")+
plotTheme()
Data source: FBI 2019, https://ucr.fbi.gov/crime-in-the-u.s/2019/crime-in-the-u.s.-2019/tables/table-8/
There are several possible reasons for the high crime rate in San Francisco:
Historical reason
Property crime rate in San Francisco hasn’t seen clear improvement in the past 10 years (2010~2020), which implies that the high crime rate is a long-term historical problem for SF.
Low arrest rate
San Francisco ranks high in crime rate while it ranks low in arrests.2 Especially the defunding of the police department causes the labour shortage, so the arrests rate cannot keep pace with the high crime rate.
Active crime
The active crimes may also account for the high crime rate. There may be more thieves in number, more active thieves or both.
“When you see an activity like that go up that much, you are either talking more active burglars, or more burglars” – UC Berkeley School of Law professor and criminologist Franklin Zimring.
Thus, it’s certainly necessary to model geospatial risk and bring predictive policing into life. In this project, a geospatial risk model is built to “borrow” the experience from the past theft crime data and forecast latent risk of thefts. There is an important underlying hypothesis that crime risk is a function of exposure to a series of geospatial risk and protective factors, such as graffiti or abandoned buildings.
A map of your outcome of interest in point form, with some description of what, when, and why you think selection bias may be an issue.
A map of theft crime is plotted with the data from DataSF. Unlike fire or property deals, theft crimes may have a much higher selection bias, which will probably leads to the systematic racism of model results.
The following conditions may result in selection bias:
If selective enforcement of thefts is the result of an officer’s preconceived or biased beliefs, then this “selection bias” will be contained into the crime outcome.
If the victim chooses not to report the crime because of modest theft value, hard or cannot reach the police, then this “selection bias” will be in the outcome.
crime.2018 = st_read("./data/crime.2018.precleaned.geojson")%>%
st_transform(crs = crs)
crime.2019 = st_read("./data/crime.2019.precleaned.geojson")%>%
st_transform(crs = crs)
nbh = st_read("./data/nbh-sf.geojson")%>%
dplyr::select(name)%>%
st_transform(crs=crs)%>%
rename(nbh.name=name)
crime.2018%>%
ggplot()+
geom_sf(size=0.3,fill="#dddddd",data=nbh)+
geom_sf(size=0.3,alpha=0.8,shape=1,color=p5[3])+
labs(title="Thefts Point Map",
subtitle="San Francisco, 2018")+
mapTheme()
As is shown in the figure above, thefts are concentrated in northwest San Francisco, where the downtown locates, and they are more sparse in the rest of the city. In addition, the nature of property thefts determine that the spatial distribution of these points is in grids, corresponding to the city road image.
Data source: DataSF, https://data.sfgov.org/Public-Safety/Police-Department-Incident-Reports-2018-to-Present/wg3w-h783
# city boundary
boundary = st_union(nbh)%>%st_buffer(1.5)
# make fishnet with cellsize 150m (around 500ft)
fishnet =
st_make_grid(boundary,
cellsize = 150, # meters
square = TRUE) %>%
.[boundary] %>%
st_sf() %>%
mutate(uniqueID = rownames(.))
# add a value of 1 to each crime, sum them with aggregate
crime_net =
dplyr::select(crime.2018) %>%
mutate(count = 1) %>%
aggregate(., fishnet, sum) %>%
mutate(count = replace_na(count, 0),
uniqueID = rownames(.),
cvID = sample(round(nrow(fishnet) / 57), # to make 100 groups
size=nrow(fishnet), replace = TRUE))
ggplot() +
geom_sf(data = crime_net, aes(fill = count), color = NA) +
scale_fill_viridis(option="C") +
labs(title = "Count of Thefts for the fishnet",
subtitle="San Francisco, 2018") +
mapTheme()
Instead of 500 meters, a 150-meter cell size is used to create the fishnet because of the smaller block size in San Francisco and denser crimes. Another reason of shrinking the cell size is that the city size is also relatively small.
Clustered spikes of high-theft-count can be observed, implying the underlying spatial process.
Point Level: Aggregate
# 311 data
var.311 = c("Graffiti","Abandoned Vehicle","Damaged Property","Homeless Concerns","Streetlights","Street and Sidewalk Cleaning")
data.311 = st_read("./data/311.cleaned.geojson")%>%
st_transform(crs)%>%
filter(service_name%in%var.311)%>%
mutate(Legend=service_name)%>%
dplyr::select(Legend)
predictors.net.1 = data.311%>%
st_join(., fishnet, join=st_within) %>%
st_drop_geometry() %>%
group_by(uniqueID,Legend) %>%
summarize(count = n()) %>%
full_join(fishnet, by = "uniqueID") %>%
spread(Legend, count, fill=0) %>%
st_sf() %>%
dplyr::select(-`<NA>`) %>%
na.omit() %>%
ungroup()
“Graffiti”,“Abandoned Vehicle”,“Damaged Property”,“Homeless Concerns”,“Streetlights”,“Street and Sidewalk Cleaning” are chosen as exposure factors according to the broken window theory.
Data Source: 311
Point Level: NN distance
nn = function(df,k=3){
return(nn_function(st_coordinates(st_centroid(predictors.net.1)),
st_coordinates(df), k=k))
}
# San Francisco Municipal Transportation stations
muniStops=read.csv("./data/muni.stops.csv")%>%
st_as_sf(wkt="shape",crs=4326)%>%
st_transform(crs)%>%
mutate(Legend="muniStops")%>%
dplyr::select(Legend)
# Bay Area Rapid Transit stations
bartStops=st_read("./data/bart.station.geojson")%>%
st_transform(crs)%>%
mutate(Legend="bartStops")%>%
dplyr::select(Legend)
# get 5 nearest neighbor
predictors.net.2 = predictors.net.1 %>%
mutate(muniStops.nn=nn(muniStops),
bartStops.nn=nn(bartStops))
# distance to downtown
downtownPoint = filter(nbh, nbh.name == "Downtown / Union Square") %>%
st_centroid()
predictors.net.2$downtownDistance =
st_distance(st_centroid(predictors.net.2),
downtownPoint) %>%
as.numeric()
# calculate nn for 311 vars
for (n in 1:(var.311%>%length())) {
vname=var.311[n]%>%paste0(".nn")
predictors.net.2 = predictors.net.2%>%
mutate(!!vname := !! data.311%>%filter(Legend==var.311[n])%>%nn())
}
Public transportation like BART(Bay Area Rapid Transit) stops and MUNI(San Francisco Municipal Transportation) stations are taken into account. All predictors are transformed into the NN(Nearest Neighbor) distance in order to improve the overall quality of the model.
Data Source:
Polygons Level: Catagories
# police districts
police = st_read("./data/police.geojson")%>%
st_transform(crs)%>%
dplyr::select(company)%>%
rename(police.ID=company)
# generaal zoning district (residential, commercial ...)
zoning.general = st_read("./data/zoning.geojson")%>%
st_transform(crs)%>%
group_by(gen)%>%
summarise()%>%
rename(zoning=gen)
# get acs data
acs.var = c(
"B02001_001E", #pop
"B02001_002E", #white alone
"B06009_005E", #bachelor
"B06012_002E", #poverty
"B19013_001E", #MedHHInc
"B25058_001E", #medRent
"B23025_007", #unemployment num
"B23025_001E" #over16
)
acs18 =
get_acs(geography = "tract", variables = acs.var,
year = 2018, state=06, county=075, geometry=T,output="wide",cache_table = T) %>%
st_transform(crs) %>%
rename(TotalPop = B02001_001E, Whites = B02001_002E,
Bachelors = B06009_005E,TotalPoverty = B06012_002E,
MedHHInc = B19013_001E, MedRent = B25058_001E,
unemployNum = B23025_007E,o16=B23025_001E) %>%
mutate(pctWhite = ifelse(TotalPop > 0, Whites / TotalPop,0),
pctBachelors = ifelse(TotalPop > 0, (Bachelors / TotalPop),0),
pctPoverty = ifelse(TotalPop > 0, TotalPoverty / TotalPop, 0),
Majority_White = ifelse(pctWhite>median(pctWhite,na.rm=T),1,0),
popDensity = TotalPop/st_area(geometry)%>%as.numeric(),
unemploymentRate = unemployNum/o16,
MedHHInc = replace_na(MedHHInc,0))%>%
dplyr::select(popDensity,pctWhite,pctBachelors,pctPoverty,Majority_White,
MedRent,MedHHInc,unemploymentRate)%>%filter(popDensity>0)
# adding Polygons
predictors.net.3 =
st_centroid(predictors.net.2) %>%
#st_join(police, by = "uniqueID") %>%
st_join(zoning.general, by = "uniqueID") %>%
st_join(nbh, by = "uniqueID") %>%
st_join(acs18, by = "uniqueID") %>%
st_drop_geometry() %>%
left_join(dplyr::select(predictors.net.2, geometry, uniqueID)) %>%
st_sf()
Census tract data of 2018 such as percentage bachelor, percentage white, population density, unemployment rate are included in the model to provide urban context.
Data Source:
ACS Tracts
Visualize the risk factors by fishnet created above.
var.names = predictors.net.3%>%colnames()%>%setdiff(c("uniqueID","geometry"))
plots1=list()
for (var in var.names){
is.discrete = predictors.net.3%>%pull(var)%>%is.character()
p = ggplot(predictors.net.3)+
geom_sf(aes_(fill=as.name(var)),size=0.001,show.legend = !is.discrete)+
scale_fill_viridis(option = "C",discrete = is.discrete, name="")+
labs(title=var)+
mapTheme()
plots1[[var]]=p
}
grid.arrange(grobs = plots1,top=textGrob("Risk Factors by Fishnet", gp=gpar(fontsize=28)))
# make polygon to neighborhoods...
crime.net.nb <- poly2nb(as_Spatial(crime_net), queen=TRUE)
# make and neighborhoods to list of weigths
crime.net.weights <- nb2listw(crime.net.nb, style="W", zero.policy=TRUE)
# local Morans'I
local_morans <- localmoran(crime_net$count, crime.net.weights, zero.policy=TRUE) %>%
as.data.frame()
# join local Moran's I results to fishnet
crime.net.localMorans <-
cbind(local_morans, as.data.frame(crime_net)) %>%
st_sf() %>%
dplyr::select(Thefts_Count = count,
Local_Morans_I = Ii,
P_Value = `Pr(z != E(Ii))`) %>%
mutate(Significant_Hotspots = ifelse(P_Value <= 0.001, 1, 0)) %>%
gather(Variable, Value, -geometry)
vars.moran = unique(crime.net.localMorans$Variable)
plots2 = list()
for(i in vars.moran){
plots2[[i]] = ggplot() +
geom_sf(data = filter(crime.net.localMorans, Variable == i),
aes(fill = Value), colour=NA) +
scale_fill_viridis(name="",option="C") +
labs(title=i) +
mapTheme(title_size = 14) + theme(legend.position="bottom")}
grid.arrange(grobs=plots2, ncol = 4,
top=textGrob("Local Morans I statistics", gp=gpar(fontsize=18)))
Figure above shows that relatively high values of Moran’s I, which represents strong and statistically significant evidence of spatial clustering of dependent variable. P-value and significant hot spot maps further shows the hot spot areas at the downtown area of SF.
final_net = left_join(crime_net, st_drop_geometry(predictors.net.3), by="uniqueID")
final_net = final_net %>%
mutate(hotspot.isSig = ifelse(local_morans[,5] <= 0.001, 1, 0)) %>%
mutate(hotspot.isSig.dist =
nn_function(st_coordinates(st_centroid(final_net)),
st_coordinates(st_centroid(filter(final_net, hotspot.isSig == 1))), k = 1))%>%
drop_na()
correlation.long= final_net%>%
st_drop_geometry() %>%
dplyr::select(-uniqueID, -cvID, -zoning,-nbh.name) %>%
gather(Variable, Value, -count)
correlation.cor =
correlation.long %>%
group_by(Variable) %>%
summarize(correlation = cor(Value, count, use = "complete.obs"))
ggplot(correlation.long, aes(Value, count)) +
geom_point(size = 0.1,color=p5[4]) +
geom_text(data = correlation.cor, aes(label = paste("r =", round(correlation, 2))),
x=-Inf, y=Inf, vjust = 1.5, hjust = -.1) +
geom_smooth(method = "lm", se = FALSE, colour = "black") +
facet_wrap(~Variable, nrow = 4, scales = "free") +
labs(title = "Theft count as a function of risk factors") +
plotTheme(title_size = 28)
A small multiple scatterplot with correlations is shown above. Important factors are damaged property, damaged streetlights, clean requests and homeless concerns.
ggplot(final_net) +
geom_histogram(aes(count),bins=50,fill=p5[4],color="black") +
labs(title = "Histogram of theft count",
subtitle = "San Francisco, 2018") +
plotTheme()
The histogram shows that the distribution of theft count per fishnet basically resembles Poisson distribution.
# define the variables we want
reg.ss.vars = "Abandoned Vehicle.nn,Damaged Property.nn,Graffiti.nn,Homeless Concerns.nn,Street and Sidewalk Cleaning.nn,Streetlights.nn,Street and Sidewalk Cleaning.nn,downtownDistance,popDensity,unemploymentRate,pctWhite,pctBachelors,MedRent,MedHHInc,hotspot.isSig,hotspot.isSig.dist"%>%str_split(",")%>%.[[1]]
reg.vars = "Abandoned Vehicle.nn,Damaged Property.nn,Graffiti.nn,Homeless Concerns.nn,Street and Sidewalk Cleaning.nn,Streetlights.nn,Street and Sidewalk Cleaning.nn,downtownDistance,popDensity,unemploymentRate,pctWhite,pctBachelors,MedRent,MedHHInc"%>%str_split(",")%>%.[[1]]
crossvalidation_id = "nbh.name"
# RUN REGRESSIONS
reg.cv = crossValidate(
dataset = final_net,
id = "cvID",
dependentVariable = "count",
indVariables = reg.vars) %>%
dplyr::select(cvID , count, Prediction, geometry)
reg.ss.cv = crossValidate(
dataset = final_net,
id = "cvID",
dependentVariable = "count",
indVariables = reg.ss.vars) %>%
dplyr::select(cvID , count, Prediction, geometry)
reg.spatialCV = crossValidate(
dataset = final_net,
id = crossvalidation_id,
dependentVariable = "count",
indVariables = reg.vars) %>%
dplyr::select(cvID=as.name(crossvalidation_id) , count, Prediction, geometry)
reg.ss.spatialCV = crossValidate(
dataset = final_net,
id = crossvalidation_id,
dependentVariable = "count",
indVariables = reg.ss.vars) %>%
dplyr::select(cvID=as.name(crossvalidation_id) , count, Prediction, geometry)
In the regression part, 4 different cross validation are used to compare the impact of factors on the final model. Specifically, two models are using spatial cross validation, i.e., to group the fishnet into neighborhoods, for each group, train the model on the other groups and test on selected model. The other two models are using groups randomly generated (about 100 groups). Besides, the variables in the models are different, spatial variables are controlled over the 4 cross validations.
# calculate errors by regression
reg.summary = rbind(
mutate(reg.cv, Error = Prediction - count,
Regression = "Random k-fold CV: Just Risk Factors"),
mutate(reg.ss.cv, Error = Prediction - count,
Regression = "Random k-fold CV: Spatial Process"),
mutate(reg.spatialCV, Error = Prediction - count,
Regression = "Spatial LOGO-CV: Just Risk Factors"),
mutate(reg.ss.spatialCV, Error = Prediction - count,
Regression = "Spatial LOGO-CV: Spatial Process")) %>%
st_sf()
error_by_reg_and_fold =
reg.summary%>%
group_by(Regression,cvID) %>%
summarize(Mean_Error = mean(Prediction - count, na.rm = T),
MAE = abs(Mean_Error)) %>%
ungroup()
## plot histogram of OOF (out of fold) errors
error_by_reg_and_fold %>%
ggplot(aes(MAE)) +
geom_histogram(bins = 30, colour="black", fill = p5[4]) +
#scale_x_continuous(breaks = seq(0, 11, by = 1)) +
#xlim(c(0,11))+
facet_wrap(~Regression) +
labs(title="Distribution of MAE", subtitle = "LOGO-CV",
x="Mean Absolute Error", y="Count")+
plotTheme()
As can be seen from the figure, there is little difference between the four CV methods. Among them, MAE of K-fold is mainly concentrated between 0 and 3 with fewer outliers, while MAE of LOGO has a higher frequency of 0, but there are more outliers, indicating that the prediction deviation of Spatial LOGO is very large for some areas, but smaller for others. The accuracy of LOGO is more extreme.
error_by_reg_and_fold %>%
filter(str_detect(Regression, "LOGO")) %>%
ggplot() +
geom_sf(data=fishnet,fill=p5[1],size=0) +
geom_sf(aes(fill = MAE)) +
facet_wrap(~Regression) +
scale_fill_viridis(option="C") +
labs(title = "Theft errors by LOGO-CV Regression") +
mapTheme() + theme(legend.position="bottom")
Just as the previous analysis, both models have large MAE in the downtown area, indicating that the models still have unexplained spatial processes. However, the spatial process cannot be explained by simply adding more spatial variables, and how to eliminate excessive MAE in downtown area is still a problem to be studied.
k = st_drop_geometry(error_by_reg_and_fold) %>%
group_by(Regression) %>%
summarize(Mean_MAE = round(mean(MAE), 2),
SD_MAE = round(sd(MAE), 2)) %>%
kable() %>%
kable_styling("striped", full_width = F) %>%
row_spec(2, color = "black", background = "#F89441aa") %>%
row_spec(4, color = "black", background = "#F89441aa")
k
| Regression | Mean_MAE | SD_MAE |
|---|---|---|
| Random k-fold CV: Just Risk Factors | 0.54 | 0.51 |
| Random k-fold CV: Spatial Process | 0.53 | 0.49 |
| Spatial LOGO-CV: Just Risk Factors | 1.21 | 2.84 |
| Spatial LOGO-CV: Spatial Process | 1.09 | 2.59 |
In general, random K-fold regression is better than LOGO regression, and spatial process plays a limited role in the model.
t = reg.summary %>%
filter(str_detect(Regression, "LOGO")) %>%
st_centroid() %>%
st_join(acs18) %>%
na.omit() %>%
st_drop_geometry() %>%
group_by(Regression, Majority_White) %>%
summarize(mean.Error = mean(Error, na.rm = T)) %>%
spread(Majority_White, mean.Error) %>%
rename(Majority_non_White=`0`,Majority_White=`1`)%>%
kable(caption = "Mean Error by tract racial context") %>%
kable_styling("striped", full_width = F)
t
| Regression | Majority_non_White | Majority_White |
|---|---|---|
| Spatial LOGO-CV: Just Risk Factors | -0.0731254 | 0.0611867 |
| Spatial LOGO-CV: Spatial Process | -0.0473177 | 0.0392849 |
In general, the prediction result of the model is larger than the actual value in the region of Majority Non White, and smaller than the actual value in the region of Majority White. Compared with common models, spatial model has smaller bias by race and better generalizability.
ppp = as.ppp(st_coordinates(crime.2018), W = st_bbox(final_net))
KD = density.ppp(ppp, 1000)
KDE_sf = as.data.frame(KD) %>%
st_as_sf(coords = c("x", "y"), crs = st_crs(final_net)) %>%
aggregate(., final_net, mean) %>%
mutate(label = "Kernel Density",
Risk_Category = ntile(value, 100),
Risk_Category = case_when(
Risk_Category >= 90 ~ "90% to 100%",
Risk_Category >= 70 & Risk_Category <= 89 ~ "70% to 89%",
Risk_Category >= 50 & Risk_Category <= 69 ~ "50% to 69%",
Risk_Category >= 30 & Risk_Category <= 49 ~ "30% to 49%",
Risk_Category >= 1 & Risk_Category <= 29 ~ "1% to 29%")) %>%
cbind(
aggregate(
dplyr::select(crime.2018) %>% mutate(count = 1), ., sum) %>%
mutate(count = replace_na(count, 0))) %>%
dplyr::select(label, Risk_Category, count)
risk_sf =
filter(reg.summary, Regression == "Spatial LOGO-CV: Spatial Process") %>%
mutate(label = "Risk Predictions",
Risk_Category = ntile(Prediction, 100),
Risk_Category = case_when(
Risk_Category >= 90 ~ "90% to 100%",
Risk_Category >= 70 & Risk_Category <= 89 ~ "70% to 89%",
Risk_Category >= 50 & Risk_Category <= 69 ~ "50% to 69%",
Risk_Category >= 30 & Risk_Category <= 49 ~ "30% to 49%",
Risk_Category >= 1 & Risk_Category <= 29 ~ "1% to 29%")) %>%
cbind(
aggregate(
dplyr::select(crime.2018) %>% mutate(count = 1), ., sum) %>%
mutate(count = replace_na(count, 0))) %>%
dplyr::select(label,Risk_Category, count)
rbind(KDE_sf, risk_sf) %>%
na.omit() %>%
gather(Variable, Value, -label, -Risk_Category, -geometry) %>%
ggplot() +
geom_sf(aes(fill = Risk_Category), colour = NA) +
geom_sf(data = sample_n(crime.2019, 3000), size = .3, colour = "#00000088") +
facet_wrap(~label, ) +
scale_fill_viridis(discrete = TRUE,option = "C") +
labs(title="Comparison of Kernel Density and Risk Predictions",
subtitle="2018 theft risk predictions; 2019 thefts") +
mapTheme()
Although the prediction accuracy of the two models is similar in the downtown area, the Risk Predictions accuracy is significantly higher in other areas while the Kernel Density is very vague for other areas.
rbind(KDE_sf, risk_sf) %>%
st_set_geometry(NULL) %>% na.omit() %>%
gather(Variable, Value, -label, -Risk_Category) %>%
group_by(label, Risk_Category) %>%
summarize(count = sum(Value)) %>%
ungroup() %>%
group_by(label) %>%
mutate(Rate_of_test_set_crimes = count / sum(count)) %>%
ggplot(aes(Risk_Category,Rate_of_test_set_crimes)) +
geom_bar(aes(fill=label), position="dodge", stat="identity") +
scale_fill_manual(values = p5[c(4,3)]) +
labs(title = "Risk prediction vs. Kernel density, 2018 Thefts") +
plotTheme() + theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
The model is satisfactory with the MAE of down to about 0.5, which means a high accuracy of result prediction. Moreover, the model has stronger generalization ability across different neighborhood contexts after adding spatial variables. As a result, the model may help control San Francisco’s high crime rate. However, the prediction accuracy of the model for some downtown areas is not high, and the reasons for this phenomenon and what kind of spatial process can be used to correct this error remain to be further explored.
However, although the model has relatively high promotion ability, selection bias is still inevitable. I think the model should not be put into use until we are sure whether the actual use of it will cause more serious racism.
1: Sanfrancisco.cbslocal.com. 2021. Fears Arise As Burglaries In San Francisco Soar. [online] Available at: https://sanfrancisco.cbslocal.com/2021/05/24/crime-and-punishment-fears-arise-as-burglaries-in-san-francisco-soar/ [Accessed 29 October 2021].